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Particle Acceleration and the Formation of Relativistic Outflows 
in Viscous Accretion Disks with Shocks 

Peter A. Becker, 1 Santabrata Das, 2 and Truong Le 3 
ABSTRACT 

In this Letter, we present a new self-consistent theory for the production 
of the relativistic outflows observed from radio-loud black hole candidates and 
active galaxies as a result of particle acceleration in hot, viscous accretion disks 
containing standing, centrifugally-supported isothermal shocks. This is the first 
work to obtain the structure of such disks for a relatively large value of the 
Shakura-Sunyaev viscosity parameter (a = 0.1), and to consider the implications 
of the shock for the acceleration of relativistic particles in viscous disks. In our 
approach, the hydrodynamics and the particle acceleration are coupled and the 
solutions are obtained self-consistently based on a rigorous mathematical method. 
We find that particle acceleration in the vicinity of the shock can provide enough 
energy to power the observed relativistic jet in M87. 

Subject headings: accretion, accretion disks - - hydrodynamics - - black hole 
physics — galaxies: jets 



1. INTRODUCTION 

It has recently been established that the acceleration of relativistic particles at a standing 
shock in an advection-dominated accretion flow (ADAF) can power the outflows frequently 
observed from radio-loud active galactic nuclei (AGNs) and galactic black-hole candidates (Le 
& Becker 2004, 2005, 2007). Radio-loud AGNs are thought to contain supermassive central 
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black holes surrounded by hot, two-temperature ADAFs with significantly sub-Eddington 
accretion rates. In these disks, the ion temperature Tj ~ 10 12 K greatly exceeds the electron 
temperature T e ~ 10 10 K (e.g., Narayan, Kato, & Honma 1997; Becker & Le 2003). The 
observed correlation between high radio luminosities and the presence of the outflows suggests 
that hot ADAF disks are able to efficiently accelerate the relativistic particles powering the 
jets. In fact, such disks are ideal sites for first-order Fermi acceleration at shocks because 
the gas is tenuous, and therefore a significant fraction of the accelerated particles are able 
to avoid thermalization and escape from the disk. Although the work of Le & Becker (2004, 
2005, 2007) was the first to establish a direct connection between the structure of the disk/jet 
system and a specific microphysical particle acceleration mechanism, their model was only 
applicable to fully inviscid (adiabatic) disks. In this Letter, we generalize their model to 
include viscous dissipation. 



2. DYNAMICAL MODEL 

The dynamical structure of the accretion disks considered here is based on the viscous 
inflow model studied by Narayan et. al (1997), Lu, Gu, k Yuan (1999), Gu & Lu (2001), 
and Becker & Le (2003), with the addition of an isothermal shock. In the one-dimensional, 
vertically-integrated, steady state ADAF scenario under consideration here, the accretion 
rate M and the angular momentum transport rate J are conserved, where 

M = AnrHpu , j = Ml-Q, (1) 

with p denoting the mass density, u the radial velocity (defined to be positive for inflow), 
H the disk half-thickness, i = r 2 Vl the specific angular momentum, Vt the angular velocity, 
and Q the torque. The vertical hydrostatic structure of the disk is described by the usual 
relations 

H(r) = ^, a 2 (r) = ^, (2) 

where a denotes the isothermal sound speed, and £k represents the Keplerian angular mo- 
mentum per unit mass for matter orbiting in the pseudo- Newtonian potential $, given by 
(Paczyhski & Wiita 1980) 

„ 2 , . GMr 3 , d$ 

with r s = 2GM/c 2 denoting the Schwarzschild radius for a black hole of mass M. The 
energy transport rate 
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is also conserved (except at the shock, if one is present), where U is the internal energy 
density, P is the gas pressure, and all quantities represent vertical averages. We assume that 
the ratio of specifics heats, 7 = (U + P)/U, remains constant throughout the flow. Note 
that the transport rates M, J, and E are all defined to be positive for inflow. 

In a steady state, the comoving radial acceleration rate in the frame of the accreting 
gas is expressed by 

Du __ du _ 1 dP t\ - t 2 
Dt dr p dr r 3 

and the torque Q is related to the gradient of the angular specific momentum t via (e.g., 
Frank, King, & Raine 1985) 

TT (dl 2£\ ar 2 a 2 

where v is the kinematic viscosity, computed using the standard Shakura-Sunyaev (1973) 
prescription, with constant a. Away from the shock location, the variation of the inter- 
nal energy density is governed by viscous dissipation and adiabatic compression, and the 
comoving rate of change of U is therefore given by (e.g., Becker & Le 2003) 

DU _ u dU _ U^dp^pu/di 2A 2 
Dt dr p dr r 2 \dr r J 

By combining various relations, one can obtain the differential dynamical equation (e.g., 
Narayan et. al 1997) 

u 2 2 7 ^ d\nu _ i 2 -t\ 2 7 /3 _ d\n£ K \ + fj-l\ u£ K (£-j) 



a 2 7 + I/ dr a 2 r 3 7 + I \r dr J \7 + 1/ aaV 

where j = J/M. This expression is supplemented by the differential conservation equation 
for the specific angular momentum, 

M = 2i _ u£ K (i-j) 
dr r ar 2 a 2 

obtained by utilizing equations (pQ) and (jSJ). Dynamical solutions are computed by simulta- 
neously integrating equations (|HJ) and ©. In order to ensure the stability of the calculations, 
the integrations are performed in the outward direction, starting with initial values near the 
event horizon computed using the boundary conditions derived by Becker & Le (2003). 

The resulting disk/ shock model depends on several parameters, namely the energy trans- 
port rate per unit mass, e = E/M, the angular momentum transport rate, j = J/M, the 
ratio of specific heats, 7, and the viscosity parameter, a. The value of e is constant in ADAF 
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disks except at the shock location. When a shock is present, we use the subscripts "-" and 
"+" to refer to quantities measured just upstream and just downstream from the shock, 
respectively. Critical points occur where the left- and right-hand sides of equation (jSj) vanish 
simultaneously. The flow must pass through at least one critical point before crossing the 
event horizon since general relativity requires supersonic inflow at the horizon. If the flow is 
smooth (shock-free), then the gas passes through only one critical point, located at radius 
r = r c . If a shock is present in the flow, then the gas passes through one critical point in 
the pre-shock region at r = r° ut , and through another in the post-shock region at r = r™ 
(Abramowicz & Chakrabarti 1990). 

The isothermal shock radius, r*, must be determined self-consistently by satisfying the 
velocity and energy jump conditions (Chakrabarti 1989) 

Ae = e + - e_ = , (10) 



w_ M 



2 ' 



where A4_ = w_/a_ is the upstream Mach number at the shock location. Note that the 
velocity jump condition we employ is slightly different from the one utilized by Le & Becker 
(2005) because those authors defined the Mach number in terms of the adiabatic sound 
speed rather than the isothermal sound speed used here. Due to the escape of energy at the 
isothermal shock location, the energy transport rate e drops from the upstream value e_ to 
the downstream value e + , and consequently Ae < 0. The power lost from the disk at the 
isothermal shock is related to the observed jet kinetic luminosity, Lj et , via 

L jet = -MAe, (11) 

where M is the accretion rate computed using the observed total energy output for a specific 
source, and the negative sign appears because Ae < 0. 

For given observational values of M, M, and Lj et , the process of determining the struc- 
ture for a disk containing an isothermal shock begins with the selection of provisional values 
for the energy inflow rate e_ and the angular momentum inflow rate j. Next we numerically 
integrate equations (|SJ) and starting from a location near the horizon and working out- 
ward towards the inner critical point. Once a solution is established that passes smoothly 
through the inner critical point, the next step is to determine the shock location r* by ensur- 
ing that the shock jump conditions (eqs. [ID]) are satisfied. However, if the resulting value 
of Ae is not consistent with equation (ITTj) . then the values of e_ and j are adjusted and 
the procedure is repeated starting with the integration from the horizon. When this step 
is successfully completed, the integration is continued from the upstream side of the shock 
outward toward the outer critical point. If the flow does not pass smoothly through the 
outer critical point, then the values of e_ and j are modified and the procedure is repeated 



- 5 - 



starting from the horizon. The end result is a unique set of values for e__ and j, along with 
the associated global solution for the disk/shock structure. 



3. PARTICLE TRANSPORT EQUATION 

In a steady-state situation, the Green's function, f G (E,r), representing the particle 
distribution resulting from the continual injection of monoenergetic seed particles with energy 
Eq from a source located at the shock radius r* satisfies the vertically-integrated transport 
equation (Le & Becker 2005, 2007) 

rr dfa 1 d / rr ^ n 9 fa 1 9 ( rr 9 f G 

-Hu-^- = —irHu) E -f^- + -— ( tHk g 



dr 3r dr dE r dr \ dr 

+ N oH E E)Hr-r,) _ r 
(AixEoYr* 

where if* denotes the disk half-thickness at the shock location, k is the spatial diffusion 
coefficient, Nq denotes the particle injection rate, c is the speed of light, and Aq is a parameter 
that describes the rate of particle escape through the upper and lower surfaces of the disk. 
The left-hand side of equation (fT2"|) represents the co-moving (advective) time derivative and 
the terms on the right-hand side describe first-order Fermi acceleration, spatial diffusion, 
the particle source, and the escape of particles from the disk, respectively. The relativistic 
particle number and energy densities associated with the Green's function are given by n(r) = 
J °° 4nE 2 f G dE and U (r) = J °° 4ttE 3 f G dE, respectively. Although the Fermi acceleration 
of the particles is concentrated at the shock, the rest of the disk also contributes to the 
particle acceleration because of the general convergence of the MHD waves in the accretion 
flow. 

After the velocity profile has been determined, we can compute the Green's function 
describing the energy and space distribution of the accelerated relativistic particles inside 
the disk by solving equation ffT2"|) . We set the injection energy using E Q = 0.002 ergs, which 
corresponds to an injected Lorentz factor Tq = Eo/(m p c 2 ) ~ 1.3, where m p is the proton 
mass. The speed of the injected particles, Vo = c(l — T^ 2 ) 1 / 2 , is about three to four times 
higher than the mean ion thermal velocity at the shock location. The seed particles are 
picked up from the high-energy tail of the Maxwellian distribution for the thermal ions. The 
particle injection rate N is computed using the energy conservation condition (cf. eq. [TT]) 

A> £o = -MAe = L jet . (13) 

This self-consistency relation ensures that the energy injection rate for the seed particles is 
equal to the energy loss rate for the background gas at the isothermal shock location. 
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To close the system, we specify the radial variation of the spatial diffusion coefficient k 
by following Le & Becker (2004, 2005, 2007), who adopted the general form 

K(r) = Kou(r)r s ^--lj , (14) 

where k,q is a dimensionless positive constant. The value of kq is determined self-consistently 
using the energy conservation relation L- ]et = L esc , where the power in the escaping particles, 
L esc , is computed using L esc = N esc E*, with N esc denoting the escape rate for the relativistic 
particles diffusing out of the disk at the shock radius r* with mean energy E* = C/(r*)/n(r*). 
The spatial diffusion coefficient k is related to the magnetic mean free path via the usual 
expression k = cA mag /3, where A mag denotes the coherence length of the magnetic field. 
Analysis of the three-dimensional random walk of the escaping particles then yields for 
the escape parameter A the result (Le & Becker 2005) A = (3k*) 2 /(cif*) 2 , where k* = 
(«_ + «+)/2. 

For values of the particle energy E > E , the source term in equation (|12p vanishes and 
the remaining equation is separable in energy and space using the functions 

■ E x v " 



L(E,r)=i — \ <p n (r), (15) 

where A n are the eigenvalues, and the spatial eigenfunctions <f n ( r ) satisfy the second-order 
ordinary differential equation 

- Hu-^ = ^-^-{rHu)(p n + ( r Hk \ _ A cH*5(r - r*) (p n ■ (16) 
dr 6r dr r dr \ dr J 

The eigenvalues A n are determined by applying suitable boundary conditions to the spatial 
eigenfunctions, as discussed by Le & Becker (2007). The eigenvalues A n and eigenfunctions 
(p n (r) can be determined numerically by computing H(r) and n(r) using equations (J2J) and 
( I14p . respectively, once a numerical solution for the inflow speed u(r) has been obtained (Le 
& Becker 2005). 

We can establish several useful general properties of the eigenfunctions by rewriting 
equation (fl6l) in the Sturm-Liouville form 



d_ 

dr 



dr 



l\ (\ n A o cS ( r ) <Pn(r) 5(r - r*) 

+ \ n uj(r)ip n (r) = — , (17) 

K[r) 



where 



,s_Su d\n(rHu) _ _ rHu 

3k dr r*H*K* 



r r 
1 s 1 s 



Ko(r*-r s ) K (r-r s ) 
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Using standard manipulations along with the asymptotic behaviors of S(r) and <^n(r)) one 
can show that the spatial eigenfunctions form an orthogonal set. It is therefore possible to 
express the Green's function / G using an expansion of the form 

f G (E,r) = J2C n ( — ) Mr), (19) 

n=0 ^ ' 

where the expansion coefficients C n are easily computed based on the orthogonality of the 
eigenfunctions <f n ( r )- 



4. APPLICATION TO M87 

In our numerical applications to M87, we adopt the observational values M ~ 3 x 10 9 M Q 
(e.g., Ford et al. 1994), M ~ 1.3 x 10" 1 M yr" 1 (Reynolds et al. 1996), and L jet = 
5.5 x 10 43 erg s -1 (Reynolds et al. 1996; Bicknell & Begelman 1996; Owen, Eilek, & Kassim 
2000). We utilize natural gravitational units, with GM = c = 1 and r s = 2, except 
as noted, and we set 7 = 1.5 (see Narayan et. al 1997). In principle, our model can 
accommodate any value for the viscosity parameter a. We set a = 0.1 here in order to 
demonstrate that shocks can exist in ADAF disks even in the presence of substantial viscosity. 
The remaining parameter values implied by the observations of M87 are e_ = 0.001516, 
e+ = -0.005746, j = 2.6257, k = 0.01632, N = 2.75 x lO^s" 1 , A> sc = 5.81 x 10 46 s _1 , 
A = 0.0153, n(r.) = 9.52 x 10 43 cm~ 3 , U(r*) = 8.94 x 10 41 ergcm" 3 , E*/E Q = 4.70, 
r* = 26.329, = 6.462, r° ut = 96.798, M- = 1.43, and = 12.10. The results obtained 
for the inflow speed u and the isothermal sound speed a in a shocked disk are plotted in 
Figure la. The value of I merges with the Keplerian value £k (eq. [3]) at r = 4, 658, which is 
the outer edge of the ADAF region. Figure la also includes the dynamical results obtained 
for u and a in a smooth (shock-free) disk with e_ = e + = 0.001516, j = 2.3988, k,q = 0.01632, 
A> = 2.75xl0 46 s" 1 , A> sc = 0, A = 0, n(r*) = 9.58x 10 43 cm" 3 , U(r m ) = 2.92xl0 41 erg cm" 3 , 
E*/E = 1.53, r c = 7.572, and = 13.61. For the purposes of comparison with the shocked 
case, we leave the source located at r = 26.329. The value of j for the smooth solution is 
selected so that £ = £k at the same radius as in the shocked disk. Our results represent the 
first dynamical solutions for ADAF disks with isothermal shocks and a significant level of 
viscosity (a = 0.1). Gu & Lu (2001, 2004) obtained solutions for ADAF disks containing 
Rankine-Hugoniot shocks, but such shocks have conserved energy transport rates, and are 
therefore not relevant for producing outflows. Note that the mean energy of the relativistic 
particles at the source radius, E*/E , is significantly larger in the shocked disk than in the 
smooth flow. 

In Figure 16 we plot the Green's function energy distribution for the accelerated particles 
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Fig. 1. — Self-consistent results for the disk structure and the particle transport obtained 
using the M87 parameters, (a) Inflow speed u(r) (solid lines) and sound speed a(r) multiplied 
by [27/(7 + 1)] 1 / 2 (dashed lines) plotted as functions of radius for shocked and smooth (shock- 
free) disks, (b) The Green's function particle distribution (eq. [19] ) plotted as functions of 
energy at various radii inside the disk in cgs units, (c) The number and energy distributions 
for the relativistic particles escaping from the disk at the shock location (eq. [20] ) plotted as 
a function of energy in cgs units. In panels (a) and (b), the heavy lines denote the shocked 
disk solutions and the thin lines represent the shock-free solutions. 
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measured at various radii inside the shocked and smooth disks. As expected, the presence of 
the shock results in a much flatter (i.e. harder) energy spectrum for the accelerated particles. 
In Figure lc, we plot the number distribution iVJ? and the energy distribution EN^ C for 
the particles escaping from the shocked disk, computed using (Le & Becker 2007) 

N e E sc (E) = (47rE) 2 r*H*cA f G (E,r*) . (20) 

The total power in the escaping particles, computed using L csc = EN^ C dE, is found to 
equal Lj Ct = 5.5 x 10 43 erg s _1 , which is a useful self-consistency check on the model. 

5. CONCLUSIONS 

In this Letter, we have obtained for the first time dynamical solutions for viscous ADAF 
disks containing shocks based on a relatively large value for the standard Shakura-Sunyaev 
viscosity parameter, a = 0.1. Utilizing a rigorous mathematical approach, we have computed 
the Green's function energy distribution for the relativistic particles accelerated in the disk, 
for both shocked and smooth disks. The results confirm that viscous disks with shocks are 
much more efficient particle accelerators than smooth disks. We conclude that the presence 
of a shock is an essential ingredient in the formation of the observed outflows. The absence 
of strong outflows from luminous X-ray sources probably reflects the fact that the gas is too 
dense to allow efficient Fermi acceleration of a relativistic particle population in the disk. This 
helps to explain the observed anticorrelation between X-ray luminosity and radio/outflow 
strength (e.g., Reynolds et al. 1996). Our results establish that the luminosity of the M87 jet 
can be understood within the context of the disk/shock/outflow model. However, collimation 
effects and radiative losses must also be considered in order to understand the subsequent 
propagation of the jet through the extragalactic environment. 
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